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A cellular automaton model is presented for random walkers with biologically motivated interactions 
favoring local alignment and leading to collective motion or swarming behavior. The degree of 
alignment is controlled by a sensitivity parameter, and a dynamical phase transition exhibiting 
spontaneous breaking of rotational symmetry occurs at a critical parameter value. The model is 
analyzed using nonequilibrium mean field theory: Dispersion relations for the critical modes are 
derived, and a phase diagram is constructed. Mean field predictions for the two critical exponents 
describing the phase transition as a function of sensitivity and density are obtained analytically. 
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When in the course of evolutionary events it became 
possible for cells to actively crawl and move towards more 
favorable habitats, this led to an acceleration of evo- 
lutionary change. Another important step was the de- 
velopment of social behavior, manifested in cooperative 
motion of individual cells or organisms. In particular, a 
change from independent crawling to cooperative motion 
is typical of life cycles in many microorganisms. 

Moving cells can orient themselves by means of indirect 
physico-chemical signals like electrical fluxes or molecular 
concentration gradients; the response of individual cells 
to such environmental information may result in collec- 
tive streaming behavior and swarm patterns. Many mod- 
els have been formulated along these lines of argument, 
all focusing on similar aspects of physico-chemical com- 
munication (see examples in Ref. Q). 

Here we are interested in the implications of direct com- 
munication between biological units (e.g. cells or birds). 
Based on the assumption that the units have an inherent 
direction of motion, and try to locally align with other 
units, several microscopic models for swarming behavior 
have recently been proposed ||-||] . These models can be 
viewed as itinerant Ay-models that can be analyzed us- 
ing renormalization group methods, starting from a pos- 
tulated equation of motion ||. 

In this Letter we take a different approach. We define a 
cellular automaton model |l that has the necessary fea- 
tures to produce swarming behavior, while the discrete- 
ness in time and space allows for relatively easy analy- 
sis. We analyze our model directly, using an approximate 
mean-field kinetic equation, and identify and derive dis- 
persion relations for the various collective modes. An 
important question is how alignment is achieved, start- 
ing from a random spatial distribution. We show that 
swarm formation is associated with a continuous dynam- 
ical phase transition, occurring when a sensitivity pa- 
rameter reaches a critical value. Spontaneous symmetry 
breaking leads to states with a global particle drift. The 
initial formation of patches is related to the fact that 



only at sufficiently large wavenumbers the density and 
longitudinal momentum modes merge to form a pair of 
propagating sound modes. We calculate the critical expo- 
nents governing the behavior of the average drift velocity 
close to criticality. 

The model we use is a lattice gas cellular automaton 
§ defined on a two-dimensional Lx L square lattice with 
periodic boundary conditions. Each node r can contain 
up to four particles in different velocity channels corre- 
sponding to nearest neighbor vectors Cj = (cos fa, sin fa) 
with fa = ir(i — l)/2 and 1 < i < 4. The state of the 
entire lattice at time t is specified by the occupation num- 
bers Si(r, t) =0,1 denoting the absence resp. presence of 
a particle in the channel (r, Cj). The state of node r is 
denoted by s(r,t) = {s t (r, i)}i<*<4- 

The evolution from time t to time t + 1 proceeds in 
two stages: first an interaction step is performed during 
which the preinteraction state {sj(r, t)} is replaced by 
a postinteraction state {cxi(r,i)} according to stochastic 
rules that are applied to each node r independently; the 
interaction step is followed by a propagation step dur- 
ing which particles move to nearest neighbor sites in the 
direction of their velocity, i.e., Sj(r + Cj, t + 1) = <x;(r, t). 

To implement the local alignment interaction we define 

4 4 

D(r,t) =^^c j; s 4 (r + c p ,i), (1) 

p=l i=l 

specifying the average flux of particles at the nearest 
neighbors of node r. We require that the number of par- 
ticles at each node, p(r,t) = p[s(r,t)] = J2t=i s.^ (r, i) , is 
conserved during interaction; this implies that the spa- 
tially averaged density of particles per node p is constant 
in time. Let J(<j) = Ci<Ti be the particle flux imme- 

diately after interaction. The transition probability from 
s(r,t) to cr(r,t) in the presence of D(r, t) is given by 

A[s <r|D] = ^S[p(a),p(s)} exp [/3D ■ 3(a)] , (2) 
where the normalization factor Z(p(s), D) is chosen such 
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FIG. 1. Swarming behavior in a cellular automaton 
model. Shown are snapshots of the systems after 100 and 
1000 time steps. Parameters are: sensitivity f3 — 1.5, system 
size L = 50, and average density p = 0.8. 



that J2 a ~~ ¥ a \^] = 1 f° r a ll s - The interaction rules 
are designed to minimize the angle between the director 
field D and the postinteraction flux J(<r). The sensitivity 
parameter (3, playing the role of an inverse temperature, 
controls the degree of local alignment: for (3 — there 
is no alignment at all; for (3 — > oo the two-dimensional 
inner product D • 3(a) — and therefore the local align- 
ment — is maximized. It will turn out that a dynamical 
phase transition occurs at a critical value (3 C of the sen- 
sitivity. Figure |l| shows the time evolution of an initially 
random distribution for [3 > (3 C . The formation of lo- 
cally aligned patches can clearly be observed. There is 
some anisotropy due to the square lattice; it is however 
straightforward to extend the model to the triangular lat- 
tice. It is an interesting question whether the phase or- 
dering kinetics shown in Fig. |l| can be described in terms 
of dynamical scaling theory H. 

To analyze the behavior of the model we consider the 
time evolution of a statistical ensemble of systems. For 
technical details we refer to Ref. B , where a model with 
only slightly different interaction rules [[l0| yet entirely 
different behavior was analyzed. In a mean-field descrip- 
tion a central role is played by the average occupation 
numbers fi(r,t) = (si(r,t)). It is assumed that at each 
time step just before interaction the probability distri- 
bution is completely factorized over channels (r, Cj), so 
that the probability to find a microstate {s.;(r)} at time 

t is given by Ur IL=i[/i(M)]' ,W [l " i*(r, i)] 1 ^ . We 
denote the factorized average by (• • -)mf- Replacing (• • •) 
by (• • -)mFj i-c, neglecting all correlations between occu- 
pation numbers, we obtain a closed evolution equation 
for /i(r,t): the nonlinear Boltzmann equation, 



/i(r + Ci,t+l) = fi(r,t) + Ii(r,t). 



(3) 



Here the term Ij(r, t) = (ai(r,t) — Si(r,t) )mf, taking val- 
ues between —1 and 1, equals the average change in the 
occupation number of channel (r, Cj) during interaction. 
It follows from the conservation of particle number, 
Ii = 0, combined with the invariance of the interac- 



tion rules under discrete rotations and translations that 
a possible solution to Eq. (||) is fi(r,t) = / = p/4. To 
assess the stability of this spatially homogeneous and sta- 
tionary solution with respect to fluctuations 5fi (r, t) = 
fi(r, t) — f we linearize Eq. (0), perform a Fourier trans- 
formation, <5/i(k, t) — J2 r e T &fi{ r it) 10' an d obtain 
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^(k,t + i)-^r J3 (k)^(k,«). 

i=i 

The mean- field or Boltzmann propagator T(k) describes 
how a small perturbation around a spatially uniform 
state evolves in time. It is given by 



r«(k) 



p=Q 



with c = and = dli(r, t)/dfj(v + c p , t)\j. It can 
be shown that 5y + fi^ = 1/4 for all this is a conse- 
quence of the fact that the outcome cr(r) of an interaction 
step only depends on s(r) through p(r) (see Eq. (||) and 
Ref. ||). For 1 < p < 4 the elements fi^ = w,j do not 
depend on p, as can be seen from the definition of D in 
Eq. (|J). We note that (ui)ij is a cyclic matrix whose first 
row has the structure (a + 7, —7, —a + 7, —7). To deter- 
mine a(/3, p) and 7(/3, p) for given values of the sensitivity 
(3 and the average density p we evaluate the expression 
(this is done numerically because of the highly nonlin- 
ear dependence on fi and /3D, combined with the large 
number of terms) 

Sj(r + ci)-/ 



(s(r+c p )} o-(r) 



/(I - /) 



x A[s a|D({s(r + cp)})] J] F(s(r + V )), 

p'=0 

where F(s) — nt=i f Si (^~ f) lsi ^ s the factorized distri- 
bution. Note that the expression for u>ij does not depend 
on r since it represents a derivative evaluated in a spa- 
tially uniform state. 

We first investigate the stability of the spatially uni- 
form state, i.e. k = 0. It can be seen that the propagator 




FIG. 2. Phase diagram for swarming model. Shown are 
the regions of stable and unstable behavior, as a function of 
sensitivity (3 and average density p. 
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ry'(k = 0) has an eigenvalue Ai = 1 with corresponding 
eigenvector e\ = (1, 1, 1, 1), reflecting the fact that the to- 
tal density is conserved. Furthermore there is a twofold 
degenerate eigenvalue \ x>y = 8a with an eigenspace 
spanned by e x = (1, 0, — 1, 0) and e y — (0, 1, 0, —1), corre- 
sponding to the x- and y components of the total particle 
flux. The remaining eigenvector e x a_ y a = (1,-1,1,-1) 
has eigenvalue \ x 2 - y 2 = I67, corresponding to the differ- 
ence between the number of horizontally and vertically 
moving particles. Numerically 7 is found to be about 
two orders of magnitude smaller than a, so that the on- 
set of instability of the homogeneous state is determined 
by the condition X x ,y = 1. The location of the criti- 
cal line in the (0, p) parameter plane is shown in Fig. ||, 
which was obtained by numerically solving the equation 
a(0,p) = 1/8. 

To see if in addition to the emergence of a global 
drift we can explain the formation of spatial structure 
in terms of the eigenvalue spectrum, we study the case 
k ^ 0. It is convenient to work with z(k) = lnA(k) so 
that excitations behave as Sf (r, t) ~ exp[z(k)£ + «k • r]. 
Unstable modes have Rez(k) > while stable modes 
have Rez(k) < 0. An imaginary part of z(k) indi- 
cates that the mode has a nonzero propagation veloc- 
ity u(k) = Im z(k)/|k|. Figure || shows that the fastest 
growth occurs at k = 0. For k ^ the degeneracy of 
X XtV is lifted, and it is then the transverse velocity (i.e., 
perpendicular to k) that grows fastest. At |k| = k p , 
with k p = fcp(k, p, 0), where k is the unit vector in 
the direction of k, the density and longitudinal velocity 
modes merge to form a pair of propagating sound-like 
modes, with Im z(k) 7^ 0, and traveling in the directions 
±k. Thus, traveling waves cannot occur on spatial scales 
larger than 27r/fc p , which may explain the length scale 
for short times of the spatial structure shown in Fig. |l|. 

Our mean-field stability analysis illuminates the nature 
of the observed phase transition. An appropriate order 
parameter is the spatially averaged velocity, 



m = ± 



which takes values between and 1. For < C we have 
p = 0. When the sensitivity parameter reaches its 
critical value, this "rest" state becomes unstable, leading 
to a breaking of rotational symmetry, and a stationary 
state where p^0. 

We have compared the results of our stability analysis 
with computer simulations. Fig. ^ shows p versus for 
averaged density p — 0.4. There is an abrupt change 
in p at ~ 0.7, which agrees well with the prediction 
0c = 0.67 obtained from our stability analysis. 

A discussion of the question whether the transition is 
first order or continuous is only meaningful if we consider 
the limit t — > 00, the analogue of the thermodynamic 
limit L — ► 00. For < C all modes are stable and we 
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FIG. 3. Eigenvalue spectrum for p — 1.6, /3 = 1.5, and 
k/x. Density (D), longitudinal (L) and transverse (T) mo- 
mentum, and sound (S±) modes are indicated. The stable 
mode that has eigenvector e x 2_ 2 at k = is not shown. 



have pit — > 00) = 0. To determine the behavior of p for 
> 0c we consider spatially homogeneous and stationary 
solutions to the nonlinear Boltzmann equation (0), i.e. 
fi(r) = fi and Ii — 0. Knowing that the "rest" solution, 
fi = f = p/4, is stable for < c (p), we expand around 
the critical point (p,0 c ): 



I i (p + Ap,0 c + A0)=J2[n 

k ^ 

+ 7 E &ikik 2 $fk 1 5fk 



da 



ik 



Oft 



A/3 ) Sf k 



k±<k 2 

where ti ikl ... kn = {d/df kl ) ■ ■ ■ (d/df k JIi @. 

We use a particular parametrization for "drift" solu- 
tions along the x-axis: 5fi — 8/3 = p, 6/2 = Sfi, and 
5fi+5f2 + Sf 3 + Sf 4 = Ap. Utilizing h = h = 0, together 
with the symmetry properties of the expansion coeffi- 
cients £li kl ... k „ and the fact that at the critical point all 
three vectors 1, c x , and c y are zero eigenvectors of f2j/-, 
we can eliminate {£/,} and for small Ap and A0 obtain 
the following equation of state: 



(cf)A0 + c p Ap)p - p? 



0. 



(4) 



Here cp and c p are positive constants that depend on the 
expansion coefficients &i kl — kn - 

Consider now the case A0 ^ 0, Ap = 0. The solution 
p = is stable only for < C , or A0 < 0. From Eq. (0) 
we see that for A0 > there is an additional, stable 
solution p ~ ^[0. Thus for the critical exponent 0' M] 
defined by p ~ (0 - 0c) 13 ' in Ref. @ we find 0' = \. A 
different exponent 6, defined by p ~ (p — p c ) s , governs 
the behavior for Ap ^ 0, A0 = 0. From Eq. (0) we 
obtain 5 = |. 
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FIG. 4. Mean velocity p, versus sensitivity (3. Obtained 
from simulation of L = 50 system at averaged density p = 1.6, 
after t — 1000 time steps. 



analyzed by observation of simulation outcomes |T(| . 
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The essential elements of our analysis are the "hydro- 
dynamic" variables density and velocity. Therefore we 
expect that many of our predictions — including the 
value of the critical exponents — should also apply to the 
continuum swarming model of Ref. ||. From a coarse- 
grained point of view our model and that of Ref. [p| are 
equivalent. In particular the noise parameter r\ in Ref. |^] 
plays a role analogous to in our model. 

Our analysis confirms the numerical finding of Ref. || 
that the phase transition is continuous, but is in con- 
flict with the results of Ref. (|]. The exponents f3' and 
S have been measured in computer simulations ||. The 
measured value (3' — 0.45 ± 0.07 is in agreement with 
our mean-field prediction f3' = |. In the case of 5 how- 
ever there is a significant deviation between the measured 
value S — 0.35 ± 0.06 and the mean-field result S = \. 
The fact that at a mean-field level (3' and 5 are equal 
supports a claim made by the authors of Ref. || that the 
observed difference between the measured values of the 
two exponents may be due to finite-size effects. 

Our model has an interesting biological interpretation 
since the dynamical phase transition suggests two possi- 
ble scenarios for a change from non-cooperative to co- 
operative behavior. On one hand, genetically caused 
minor microscopic effects on receptor properties of in- 
teracting cells influencing their sensitivity can have se- 
vere macroscopic implications with respect to swarming 
if they occur close to criticality (cf. Fig. 2). On the other 
hand, a transition from the stable into the unstable re- 
gion can also be achieved by simply increasing cell density 
(cf. Fig. 2). This result provides a possible clue to ex- 
plain the behavioral change between non-cooperative and 
cooperative stages in individual life cycles of some bac- 
teria and amoebae in which a reproductive feeding phase 
of individually moving cells is followed by social (coordi- 
nated) aggregation. Other models for complex bacterial 
pattern formation, including vortex and colony organi- 
zation, have also been proposed (see Refs. @,|l5| and 
references therein). 

Finally, we want to stress that the methods employed 
here can easily be adapted to gain theoretical insight in 
the behavior of a wide range of biologically motivated 
cellular automaton models, that so far have mainly been 
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